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Using a viscous hydrodynamic model coupled to a hadronic cascade code, numerous features of the 
dynamics and equilibrium properties are explored for their impact on femtoscopic measurements. 
The equation of state, viscous parameters and initial conditions are investigated. We find that 
femtoscopy is affected by numerous model features at the 10% level, and that by including features 
and adjusting unknown parameters, one can explain experimental source size measurements to better 
than 10%. 



I. OVERVIEW 



The bulk properties of QCD matter, as created in relativistic heavy ion collisions, largely manifest themselves 
in soft hadronic observables of particles with transverse momentum less than one GeV/c. These observables 
can be divided into three classes: spectra, flow (or large-scale correlations), and correlations at small relative 
momentum. This last class is referred to as femtoscopy [T] since these correlations are used to determine 
space-time characteristics of emitting sources. Correlation functions, C(P,q), can be linked to the outgoing 
phase-space distributions, or more precisely the source function ^(Pjr), which describes the probability that 
two particles with the same velocity, whose total momentum is P, are separated by r in their asymptotic 
trajectory. Due to their inherent six-dimensional nature, correlations have proven to be the most difficult of 
all RHIC observables to fit with full dynamic models. The measurements are amenable to being fit by simple 
geometric models of the final state, provided that the models incorporate strong radial collective flow, and a 
rapid dissolution into a thermal assortment of resonances [21 El 131 [5] . However, many dynamic models, especially 
hybrid hydrodynamic/cascade descriptions, lead to more extended emission durations which lead to signicantly 
different shapes for the outgoing phase space distributions. 

The information in correlations are often reduced to Gaussian source parameters, i?out, -Rsidc and i?iong, which 
are functions of the transverse momentum kt, and describe the shape of the outgoing phase-space distribution 
of zero-rapidity particles with a specific kt and a specific azimuthal angle. Here, i?iong refers to the longitudinal 
dimension, along the beam axis, -Rout describes the outward dimension, parallel to the momentum, and i?side 
refers to the sideward dimension, perpendicular to the beam and to the particle's velocity. Asymptotically, the 
Gaussian form fits the phase space density to the form. 
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The term Vpt -|- a is irrelevant for identical particles since correlations are only sensitive to the relative position 
of two particles of the same velocity. For non-identical particles, one would also be sensitive to the relative 
position of centroids for the two species, ai — 02. 

This study focuses on how these parameters are affected by choices made in modeling the reaction. The 
femtoscopic relation to the equation of state has long been studied. First, a stiffer equation of state leads to 
more rapid expansions, with emission at earlier times and more confined to a brief burst. The reduction in 
the mean emission time reduces i?iong as the system has less time to expand longitudinally before emission. 
The increased suddenness leads to a shorter -Rout relative to -Rsidc, as long-lived emission allows those particles 
emitted earlier to move ahead of the later-emitted particles along the outward direction. Furthermore, a softer 
equation of state leads to higher entropy, and for fixed spectra, the entropy will be grow for increasing source 
volumes, V ~ -Rout-Rside-Riong- Sincc the total entropy can be ascertained in a quasi-model-independent fashion 
from spectra and source dimensions, and since entropy is conserved during the expansion, the product of the 
three dimensions is strongly linked to the equation of state 

Femtoscopic source sizes are also affected by non-equilibrium aspects of the dynamics. Bulk viscosity, which 
is expected to be significant in the neighborhood of the critical temperature when the system struggles to 
maintain equilibrium, lowers the effective pressure and thus increases the entropy and leads to larger source 
dimensions. Shear viscosity is mainly important at early times, when velocity gradients are large and highly 
anisotropic. This leads to an enhancement in the transverse components of the energy tensor, which accelerates 
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the transverse expansion and gives smaller values of i?iong and i?out relative to i?side [TJ [S] . At the earliest times, 
before even viscous hydrodynamics is applicable, the pre-equilibrium state might be dominated by longitudinal 
color fields. These fields can, in principle, lead to exceptionally strong transverse components to the stress- 
energy tensor, which would amplify the effects of shear viscosity. The importance of early acceleration in 
explaining the experimental i?out/^sidc ratio has also been studied by incorporating initial transverse flow into 
ideal hydrodynamics and by adding a strong repulsive potential into microscopic cascade codes [10]. Free- 
streaming during the first fm/c increases the transverse pressure relative to the longitudinal pressure, which 
increases radial flow more than elliptic flow This had been thought to make it difficult to simultaneously 
fit both spectra and elliptic flow, though this was accomplished in jT^]. This issue might be resolved by better 
understanding the interface between the initial-state description and the following hydrodynamic description. 

The initial density proflle also affects acceleration at early times [T^] . Within the typical nuclear cross-section 
(40 mb), a single nucleon will interact with multiple nucleons from the opposite beam. Depending on the 
theoretical picture, e.g., color-glass-condensate-inspired or wounded nucleon, the average radius of the initial 
flreball can vary by ^ 10%, with a more compact source being more explosive and leading to smaller values of 

-^out/-^sidc- 

To investigate these effects we apply a relativistic viscous hydrodynamic model, which couples to a cascade 
code for modeling the hadronic breakup stage and the decay of outgoing resonances. The outgoing phase space 
points for pions are used to generate source functions, and after convoluting with the squared wave function, 
generate two-pion correlations. These are then treated as data and fit to correlation functions from Gaussian 
sources to extract i^out, -Rside and i?iong as a function of kt- The hydrodynamic code uses Israel-Stewart 
equations which are modified to allow one to tune to the anisotropics of the initial stress-energy tensor. Both 
the hydrodynamic and cascade descriptions are built on an assumption of azimuthal symmetry and boost- 
invariance. This prohibits a simultaneous analysis of elliptic fiow, or a study of longitudinal acceleration, which 
is known to affect results at the 5-10% level. In addition to investigating all the sensitivities alluded to above, 
we compare data from STAR collaboration at the Relativistic Heavy Ion Collider (RHIC) . We do find solutions 
which come close to the the data without employing any particularly disquieting assumptions or any parameters 
outside what we would consider a reasonable range. Although this paper focuses on femtoscopy, the mean pt of 
various calculations is also presented and compared to data. 

After reviewing details of the model in the next section, the following section is devoted to the effects of varying 
the equation of state, viscosities, and initial conditions. The summary is devoted to drawing conclusions with 
an emphasis on understanding what future improvements in models and analysis are needed to reach rigorous 
quantitative statements about the microscopic properties of the QCD matter formed in relativistic heavy ion 
collisions. 



II. THE MODEL 

For generating interferometric source functions, phase space points are first generated from a viscous hydro- 
dynamic model, then fed into a cascade model which models the low-density hadronic stage of the collision. 
Both models are written in terms of the variables r, ry, r and (j), where r = ^Jt"^ — is the proper time, 
r\ = sinh~^(z/T) represents the longitudinal position, and r and (\) represent the radial position and azimuthal 
angle. Both models were developed assuming radial symmetry and boost invariance which eliminates r\ and <^ 
from consideration. By reducing the dimensionality, both speed and accuracy are vastly improved. The viscous 
hydrodynamic model is based on the formalism in |13l . with more details provided below. The subsequent 
subsections provide details of the three model components. 

A. Viscous Hydrodynamic Model 

First, a review of the modified Israel-Stewart formalism described in [13 is presented. A basic description of 
viscosity and the Navier-Stokes equation can be found in [2]. Recently, Israel-Stewart hydrodynamics [15j has 
been extended and applied to nuclear physics [111 [13 Ell Ell HH]- In ideal hydrodynamics, the stress energy 
tensor becomes Pbij when viewed in the fluid rest frame (Here latin indices refer to spatial components only). 
Viscous hydrodynamics deals with the deviation of Tij from Pbij . For all the hydrodynamic calculations here, 
the fluid rest frame is defined such that Tqi = 0, and diffusion of conserved particle numbers through fluid 
elements is ignored. In the fluid frame the deviations of T!y can be expressed through flve independent traceless 
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components, a^, and the deviation of the trace, b. 



b = -{T,,+Tyy+T,,)-P, (2) 
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The shear components ai are related on a one-to-one basis to the five velocity gradients, LOi, 



LOl = d.^V^ - dyVy, (3) 

1 

W3 = {d^Vy + dyv.^) , LOi = {dxVz + dzVx) , ^5 = {dyVz + dzVy) . 



u^i = {d^v^ + dyVy - 2dzVz) , 



With these definitions, the Navier-Stokes equations become 

ai^-rjLOi, b=-(y-v. (4) 

For Israel-Stewart equations of motion, and 6 are not fixed as is the case for Navier-Stokes equations, but 
instead are dynamic objects. The ratios aijoa and bja^ should decay exponentially toward the Navier-Stokes 
values, where Oa and (Jb are related to the fluctuation of the stress energy tensor at fixed energy density. 



-2 _ / j3 
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d''r(6(0)6(r)), (5) 



,2 _ / j3 



d^r(a,(0)a,(r)), 

where no sum is implied in the last expression. To restrict the values of ai and b to ranges ±amax and ifemax 
respectively, the Israel-Stewart equations are modified by mapping a.^ and b to yi and x through hyperbolic 
tangents. The variables yi and x will follow the Israel-Stewart equations above and can become arbitrarily 
large, while ai and b will be restricted. 

(yi - Wi/c^a) , (6) 

Tr, 



dt 



(^ay 



y ^ \Jvi+vh 



Ui / 9 , 9 

ai = a — , a — y af + 

dx 1 , , , 

— = (x - Cv • v/at) , 

at Tf, 

b = 6max tanh ( 

As derived in [T3] and [21], the lifetimes, fluctuations and viscosities are not independent, 

C-^P. (7) 

The equations of motion are solved by storing the velocities and energy densities in a mesh defined by the 
radial coordinate. Following the ideas of the one-dimensional calculation in [22], mesh points are not stored at 
equal times but at varying times that enforce local simultaneity, i.e., u ■ Ax = 0, where Ax is the four vector 
describing the separation of two neighboring mesh points. Using the integrated distance along the mesh as seen 
by comovers, 

d£, de = ^/-[dx~u{u-dx)]^, (8) 
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the acceleration in the fluid frame, which equals the rate of change of the transverse rapidity, takes on a simple 
form, 

--y, = a{£) = - . (9) 

UT J-00 + -Ixx 

In order to maintain simultaneity between neighboring mesh points, the time step depends on £, 

dT{e) = ST{e = 0) exp 1^ d£' a{£') | . (10) 

To complete the equations of motion, an expression is needed for the evolution of the energy density. This is 
done by considering the change in the internal energy within a cell defined by adjacent mesh points and a fixed 
small rapidity range Srj. In the fluid frame, the volume of the cell is 

Ay = {T5ri)2'KRM, (11) 

where R is the radius as viewed in the laboratory frame. In a time step dr the volume increases both due the 
increase in the longitudinal dimension d{TSri) and by the increase in the transverse dimensions. Writing 

dAV = dAV^ + dAV,, (12) 
dAK = {TSr])d{2TrRAe), dAV,_ = (2'KRA£)d{Tdr]), 

the change in the internal energy of the cell is 

dAU ^ ~T,,dAV^ ~ T,,dAV„ (13) 

where z is the longitudinal direction and x refers to the radial direction. Given the internal energy and volume, 
one then knows the local energy density e which closes the equations of motion. For the ideal case, Tij = P5ij, 
one recovers dU = —PdV, which implies entropy conservation. Indeed, when the code was run in this limit, 
entropy was conserved to better than 0.2 percent. 

The equation of state used for the runs shown here consisted of three parts. For temperatures below 170 MeV, 
the equation of state was that of a hadronic gas. For a given cell, the pressure was calculated as a function 
of the energy density and the density of five conserved charges. The conserved charges were the number of 
strange plus anti-strange quarks, the number of baryons plus antibaryons, the effective number of pions (e.g., a 
p meson counts as two pions) , the number of rys and the number of los. Only the standard octet mesons and octet 
and decuplet baryons were considered, i.e., the tt, K, ri,p,uj, K* ,7]' ,uj, 4> mesons and the p, n. A, E, S, A, E*, S*, 
baryons. The details of which particle numbers were fixed was not particularly important because the breakup 
density was chosen to be 400 MeV/fm"^, which allowed the cells very little time to adjust their chemistry before 
the evolution was taken over by the cascade code. 

For an intermediate range of energy densities, e/i < e < + L, the equation of state was chosen to have a 
constant speed of sound, i.e., P = Ph + c^ixod(^ ~ ^h)- Here, e/i is the energy density of an equilibrated hadron 
gas with a temperature of 170 MeV. In the limit of Cmixcd = 0, the equation of state becomes that of a first 
order phase transition with latent heat L. For energy densities above e/j + L, the speed of sound was bumped 
up to (? = 0.3 to be consistent with lattice gauge theory ^34j. The simple form for the equation of state was 
used so that by varying L and c,„ixcd one could study the sensitivity to the equation of state. 

The ratio of the shear viscosity to entropy was fixed above T^. According to the KSS conjecture, this ratio 



should stay above l/47r [23]. Results for varying 77/s are investigated in Sec. Ill The fluctuation cr^ was 
calculated by considering fluctuations of a massless gas. This gives cr^ = (47r/5)T^s. The relaxation time 
from Eq. Q is then Tq = {5/4:Tr)r]/Ts. For e < Ch, the relaxation time was chosen as Ta — l/{na), with 
(T = 25 mb and n being the hadron density. This energy-averaged cross-section for a hadronic gas would be 
somewhat higher than 25 mb, but much of that cross section would be more forward peaked which reduces the 
effectiveness of collisions to thermalize the matter. It would not be surprising if more sophisticated calculations 
of the relaxation time would differ by a few tens of percent. The fluctuation aa was determined by considering 
the fluctuations inside a hadron gas. 



a' 



|dV T.,(0)T.,(r) = 1 Y: E (2ja + l)/^/.(p)^. (14) 

particles i * species a 



The relaxation time is then given by Eq. ([7|. For the intermediate region, e/i < e < eh + L, both 77 and cr^ were 
chosen to vary linearly with the energy density from the hadronic value at th to the value for the lower end of 
the plasma region. Relaxation times were then chosen according to Eq. Q. 
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Bulk viscosities are expected to be negligible except near Tc [21] ■ This can be understood by considering an 
isotropically expanding thermalized g 3jS , I.e., cl Hubble expansion, which for photons maintains a thermal form 
to the photon spectrum with the temperature falling ~ l/r. For a non-relativistic gas, a thermal form also 
ensues, but with the temperature falling as If the shape of the momentum-space distribution is already 

thermal, collisions are needed to maintain the thermal value for Tu. Bulk viscosities thus disappear high above 
Tc where the temperatures far exceed the quark mass, and for temperatures well below the pion mass. 

In contrast, near the degrees of freedom and the condensed fields need to change in order to maintain 
equilibrium, which leads to a bulk viscosity in that region [24, 25 . The bulk viscosity was chosen to be zero 
for e > e/i + L and for e < eh- In the middle of the mixed region, e = + L/2, ( was set to a maximum 
value, Cnax- To make the bulk viscosity a continuous function, it was chosen to linearly fall with energy 
density above and below the maximum value so that it returned to zero at the boundary of the mixed region. 
Arbitrarily, the relaxation time was chosen to be 5?i/(47rT), which equals the minimum relaxation time for the 
shear in the plasma phase if one is at the KSS limit. The treatment of the effects of bulk viscosity here are 
undoubtedly naive. Since the principal source of bulk viscosity might be the non-equilibrium chiral condensate, 
or cr field, relaxation times might be very large, the response might be very non-linear, and the behavior 
might be oscillatory, which contradicts the Israel-Stewart assumption that non-equilibrium deviations decay 
exponentially. Thus, the investigations can probably only qualitatively describe the impact of non-equilibrium 
effects towards the trace of the stress-energy tensor. Bulk effects due to non-equilibrium fields would be better 
treated by a simultaneous solution of the respective wave equations coupled to the hydrodynamic medium. 

Figures [l] and [2] display the hydrodynamic evolution for the default parameter set. The kinks in the collective 
transverse rapidity shown in Fig. [T] arise from the region near Tc- As the matter expands into a region with 
lower speed of sound a pulse builds up similar to a tsunami. Bulk viscosity, which lowers the effective pressure, 
(Tii) = (1/3) (Tea; -|- Tyy + T z z) ^ lu thls region amplifies the pulse. The pulse largely dissipates by the time final 
breakup occurs as the rapidity profile becomes linear and the energy density profile also becomes smooth. 

Viscous effects on the stress energy tensor are illustrated in the three panels of Fig. [2] The effective pressure 
{Tii) is displayed in the upper panel. Since the bulk viscosity is set to zero outside the intermediate region, 
the ratio {Tii)/P varies from unity only in this region. If not for the saturation enforced by Eq. (|6]), the 
effective pressure might fall below zero. The size of the effect is enhanced by the pulse which results in large 
velocity gradients at the boundaries of the pulse. The anisotropy of Tij is shown in the lower two panels. 
At tq the anisotropy was inserted as a boundary condition with Tzz set to zero, which equivalently gives 
{Tzz — {Tii))/P = — 1 as shown in the middle panel. This is also the saturated value as enforced by Eq. ([6]), and 
is maintained for some time due to the large velocity gradients at early times. At the edge of the fireball, where 
the matter is in the lower density hadronic phase, the strong anisotropy remains due to the large viscosity at low 
density. However, the behavior in this region is somewhat irrelevant as it is below the breakup density and is 
handled by the cascade description described in the next section. The lower panel shows T^x — Tyy, which differs 
from zero mainly near Tc due to the radial pulse. The large variations shown in Fig. |2]shows how Navier-Stokes 
treatments, which can lead to arbitrarily large deviations of the stress-energy tensor, are questionable at early 
times or in the region of the radial pulse. 

The hydrodynamic module was run until the entire system fell below the breakup density. A sampling of 
emitted hadrons was generated from the entire evolution. At each time step, particles were generated from 
thermal surface emission of the outermost cell whose energy density was above the breakup density. The 
generation of particles was consistent with the temperature, density, and the anisotropy of the stress-energy 
tensor. When a cell's energy density fell below the breakup density, particles were emitted from that cell 
in the same manner, except according to volume emission. In order to make such emission consistent with 
the anisotropy of the stress energy tensor, the particles were first generated according to an isotropic thermal 
distribution. The momenta p^, Py and pz as seen in the fluid frame were then scaled by factors A^;, Ay and Xz 
respectively, where 

A, = Vryp. (15) 

The same mechanism is used for surface emission, along with an additional factor taking into account the rate at 
which the particles leave the surface. For a given particle moving along a collision-less trajectory, the momenta 
as measured in the local rest frame should scale inversely with the time between collisions. For each component 
of the momentum, piocai = PoTq / {t + Tc) , where tq is the inverse velocity gradient at the time of the last collision 
along the given direction and Tc is the time since the last collision. For non-relativistic particles one can derive 
the simple scaling form for the various momenta shown above. However, for lighter particles the simple scaling 
is only approximately justified. Future versions of the program will apply a more sophisticated mechanism for 
generating particles. 
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FIG. 1: (color online) The transverse rapidity (upper panel) and energy density (lower panel, multiplied by r) profiles 
for three times: r = 1 fm/c (solid line), r = 3 fm/c (dashed line) and t = 6 fm/c (dotted line). The lower speed of sound 
near Tc causes a tsunami-like pulse to grow which then largely dissipates before breakup. 



B. Hadronic cascade Model 



For energy densities below 400 MeV/fm^, a cascade code is used to describe the evolution. The cascade 
simulates the evolution of the particles as straight line trajectories, punctuated by collisions whose probability 
is determined by a combination of a fixed cross section of 15 mb, along with resonant absorptions and decays. 
The resonant cross sections use a simple Breit-Wigner form with fixed lifetimes, and all collisions and decays 
axe treated as s-waves. Only resonances from the standard meson octets or from the baryon octet and decuplet 
are included. This is a simple treatment, with no mean field or Bose effects, but should provide a sufficiently 
reasonable description of the breakup stage for the interferometric studies presented here. 

Particles are entered into the cascade description from a Monte Carlo list generated by the hydrodynamic 
module. Along with the list of particles, the cascade module is also given a description of the position of the 
emitting surface as a function of time. Any particle that returns to the interior of the surface during the cascade 
description is deleted. Assuming that the hydrodynamic code, with its inclusion of shear anisotropics, accurately 
models the behavior of a hadron gas for energies near the 400 MeV/fm^, this should provide a fully consistent 
interface. For all of the parameter sets studied here, the hyper-surface from which particles are created by the 
hydrodynamic module rapidly collapses resulting in time-like emission for the vast majority of particles. The 
percentage of particles that are re-absorbed into the hyper-surface during the cascade is only about one percent. 

One potential issue with many cascade codes is that the finite interaction range, i.e., particles collide at a 
finite interaction range of approximately y/a/n, leads to viscous effects [IB] . Usually, such effects are minimized 
by oversampling the distribution by a factor iVsampie which reduces the cross sections by 1/-/Vsampie and the 
interaction ranges by 1/ -/Vsampic • However, in the cascade description applied here, the interaction range is 
set to zero, thus eliminating such numerical viscosities. This is accomplished by exploiting boost invariance and 
azimuthal invariance, which allows the trajectories to be treated as radii evolving as a function of the proper 
time, t'(t). When two particles have the same radius, a probability is calculated for their colliding given the 
cross section and the fact that the sampling covers 27r radians and one unit of rj. Given that the two radii are 
equal and that the other coordinates are irrelevant given the symmetries, the effective interaction range is zero. 
Furthermore, since any correlations from collisions or resonant decays are spread out over a wide range of r] 
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FIG. 2: (color online) Deviations of the stress-energy tensor due to viscous effects are displayed as a function of r for 
three times: r = 1 fm/c (solid line), r = 3 fm/c (dashed line) and r = 6 fm/c (dotted line). The effective pressure, 
(Tii) = {T^x + Tyy + Tzz)/3 scaled by the pressure (upper panel) deviates from unity due to bulk viscosity which is only 
non-zero near Tc- The longitudinal components (T^z — {Tii))/P begins at the saturation value of -1 enforce by Eq. ([6|. 
The ratio moves toward zero as the velocity gradient lessens, but less so for large r due to the large viscosity at low 
density. The lower panel shows that an anisotropy in {T^x —Tyy)/P grows mainly in the region where the outgoing radial 
pulse creates a pulse in the radial component of the velocity gradient. 



and (j), this treatment should come extremely close to a true Boltzmann description even though particles are 
represented on a one-to-one basis. 

The algorithm is optimized by storing the information for each particle in a list ordered by radius. A second 
ordered list stores the list of pairs that will cross and is ordered by crossing times. The crossings are executed 
in order of time, and since the list is ordered, new crossings need only be calculated for the nearest neighbors 
of those particles that have crossed. When two particles cross, one needs to calculate the probability that they 
collide, or merge to form a resonance. This is related to the ratio of the cross section to the area of the cylinder 
over which the particles are spread, 27rrT (assuming one unit of rj is being modeled). The exact probability is 
complicated by the relative angles of the particles and relativistic effects, and is given by: 



Scattering Probability = - — — — , (16) 

2ttRt Poq,. - P^go 

Pa did ol ( \ol t^ol ^ (.-f 1 ) 
= q^{pi~P2) ~P ^ • 

Here x refers to the radial components. 

All collisions before 25 fm/c are simulated, with decays being performed until they are exhausted. Weak 
decays are allowed to take place, except for charged kaons, pions and iiTiong mesons. The point at which each 
particle had its last collision is recorded to be used for calculating spectra and femtoscopic correlations. A 
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FIG. 3: (color online) Positions of last interaction for pions with = 300, Py = Q MeV/c. The outward position has a 
modestly positive correlation with emission time. This correlation is indicative of an exploding source. Even though the 
duration emission is rather long, over 10 fm/c, the correlation allows Rout/Raidc to be close to be close to unity. 



sampling of phase space points is displayed in Fig. |3]for particles with transverse momentum pt — 300 MeV/c, 
and reveals a modest positive correlation between the outward position and time. This is opposite what one 
expects for an inwardly burning source whose emitting surface would move inward with time. Instead, it suggests 
an exploding source. 

A glance at Fig. [3] reveals a strikingly different source than that extracted from Blast-wave models. First, 
both the average lifetime and the duration of the emission are longer than those determined by blast-wave 
models [21 H]- Secondly, whereas the blast- wave analysis of [2] suggested a negative x — t correlation, there is 
a modest positive correlation in Fig. |3] Such a positive correlation was also seen in the AMPT model [57], 
which is based on a cascade picture for both partons and mesons. Since the hyper-surface representing the 
transition from hydrodynamics to the Boltzmann approach has a negative x — t correlation, this emphasizes the 
importance of accurately accounting for the breakup dynamics with a microscopic model. An underestimate of 
the emission duration for blast-wave models is expected if the blast-wave model employed an a; — i correlation 
of the wrong sign, as a positive correlation allows one to maintain a small i?out/^sidc ratio despite a longer 
emission duration. The mean emission time is associated with the ratio -Riong/i^z,thorm, where Wz, therm is the 
thermal velocity for longitudinal motion. Since blast wave analyses typically ignore shear effects and thus assume 
thermal motion is locally isotropic, they would overestimate t^z, therm if in fact the local momentum distribution 
is broader in the transverse plane than along the longitudinal direction. An overestimate of Uz. therm would lead 
to an underestimate of the lifetime. 

The algorithm is both efficient and accurate. The procedure eliminates artifacts associated with particles 
interacting at a finite separation, and running on a single CPU, an event is performed in less than ten seconds. 
However, the approach has one numerical disadvantage. When calculating the value of r = \Jt'^ — at which 
two particles will reach the same radius, one must solve quartic equations. Numerical errors for such solutions 
are non-negligible which occasionally lead to particles being propagated in such a way that violates the ordering 
by radius. This calculation is performed millions of times within a single event, and the violations tend to occur 
approximately once per every ten events. In such an instance, the event is abandoned. This does not seriously 
detract from the numerical efficiency, but such failures significantly complicated the construction of the code, 
as frequent error checking is required. 
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C. Generating and Fitting Correlation Functions 

Correlations can be generated via the Koonin equation, 

C(P,q) = J d\ S{P/2,r)\<j>{q,v)\\ (17) 

c..p ^ ^ / d^Xgd^xt s{P/2, Xa)s{P/2, xt)S{r - (x; - x'J) 
^ ' ' Jd*Xgd*x,s{P/2,Xg)s{P/2,Xk) 

where P is the pair's momentum, r is the spatial separation of the particles in the frame of the pair, and (f> 
is the outgoing relative wave function. The probability of emitting a pion of momentum p from space-time 
point X is s(p,x), with x' being the coordinate in the pair frame. The source function S'(P,r) is simply the 
normalized probability that two particles of the same momentum, P/2, are separated by r in the pair's center 
of mass. The relative wave function incorporates quantum symmetrization, and both the Coulomb and strong 
interaction between the two pions. For the calculations shown here q refers to one half the relative momentum 
as measured in the pair frame. 

In practice, the Koonin equation is straight-forward to implement. To calculate C(P,q), one first extracts 
the subset of phase-space points from the output of the Boltzmann codes whose transverse momenta are within 
5 MeV/c of P/2. For every pair in the subset, one calculates \(f>{q_, r)p for an array of q values. The same set of 
pairs is used for every value of q. The correlation function is then the average of \(f>\'^ for the pairs. Statistics for 
such calculations are greatly enhanced by the rotational and boost invariances, as every particle's phase space 
points can be rotated and boosted so that it has zero longitudinal momentum and travels in a given azimuthal 
direction. 

If one neglects the inter-pair Coulomb and strong interactions, the calculations can be greatly accelerated by 
calculating 

p(q)=^e2''i-% (18) 

i 

where the sum covers the N particles in the subset used above, and is the position of the z*'' particle. For 
large N, correlations can be generated by simply squaring the sum, 

C(q) = l + ^|p(q)P. (19) 

Since one never has to evaluate a double-sum, this method is quicker than the alternate method. Although 
it can only be applied if one neglects strong and Coulomb forces, this method should be sufficient if one is 
generating correlations for the purpose of finding effective Gaussian source sizes. Both methods were tried for 
the calculation with the default parameters, with the comparison being illustrated in Fig. [4] Since the differences 
in the extracted Gaussian source sizes were small and the trends of interest are unlikely to be affected, the latter 
method was chosen. Although the calculations using the full wave functions are more realistic, it should be 
pointed out that experimental analyses have generated source radii by dividing the experimental correlation 
function by a g— dependent factor with the purpose of dividing away the effect of the Coulomb force in affecting 
correlation funcitons [281 EH] . Since the Coulomb correction factors are based on isotropic Gaussian sources, 
the procedure is not exact. Now that the discrepancy between models and experiments are less than 10%, the 
errors introduced in this procedure should be re-examined. In particular, errors should be checked for radii at 
higher kf. For kt ~ 500 MeV/c, source functions are highly anisotropic in the frame of the pair due to Lorentz 
contraction, with i?out approaching five times i?sidc, whereas the correction factors are built assuming that the 
source shape is isotropic in that frame. 

To generate Gaussian radii, correlations were calculated on a three-dimensional mesh in q. These were then 
compared to predictions for C(q) for Gaussian sources. The radii were chosen to minimize the sum of the 
squared radii. Even though the correlations were remarkably non-Gaussian for q < 10 MeV/c, the radii were 
remarkably robust, and did not change appreciably if one neglected the low q points in the fit. Since the mesh 
and q values were generated in the pair frame, the outward source size was then Lorentz contracted so that it 
represented the shape of the outgoing phase space density in the laboratory frame. 



III. FEMTOSCOPIC RAMIFICATIONS OF ADJUSTING THE EQUATION OF STATE, 

VISCOSITY AND INITIAL CONDITIONS 

One of the prime motivations for interferometric measurements was the possibility of observing a long-lived 
mixed phase, which could only happen if the equation of state was first order with a large latent heat. Some 
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FIG. 4: (color online) Gaussian source sizes were found by fitting generated correlation functions to those from Gaussian 
sources. For the default parameter set, correlations were calculated from Eq. (17 1 with (triangles) and without (circles) 
the effects of Coulomb and strong interaction in the relative wave function. Since calculations are much quicker without 
the interactions, and since differences are small, interactions are neglected for the calculations in the next section. 
The differences between the default calculation differs from experimental values (stars) (251 are significantly lessened as 
compared to calculations based on ideal hydrodynamics [301 131| . 



bag model parameterizations employed in the early days of the field had latent heats of several GeV/fm'^. If 
that were the case, extracted lifetimes at the AGS and SPS might have been several tens of fm/c depending 
on the amount of initial stopping [32l |33] . The longest lifetimes would occur for conditions where the interior 
energy density was initially at the peak value for the mixed phase. Since a mixed phase has zero sound velocity, 
there would no impetus for explosion, and instead the outside would emit hadrons like a burning log. Lattice 
calculations now preclude such equations of state, and indeed, no such long lived phases have been observed. 
Instead, in lattice calculations the speed of sound appears to dip down to about ~ 0.1, before re-stiffening 
to 0.3 at high temperatures [34^ By including resonances in a hadron gas, the speed of sound is expected 
to be ^ 0.15 below Tc- Given that initial energy densities at RHIC are well above those of the soft region, 
one should expect RHIC collisions to be more explosive and shorter-lived that those at the AGS and SPS. 
Qualitatively, these expectations have been met. However, quantitatively describing the source sizes with full 
hydrodynamic models has proved elusive. Femtoscopy provides a six-dimensional test of any dynamical model, 
so it should not be surprising that reproducing experimental source sizes requires using a realistic equation of 
state, accurately modeling viscous effects and using correct initial conditions. We explore the impact of each of 
these three aspects of the modeling in the next three subsections. Results from the default parameter set are 
compared to results where an isolated parameter set has been adjusted. For each calculation, the initial energy 
density is adjusted so that the final dNch/drj ~ 690 |55] . 

Radial flow, and in turn spectra, are also affected by all the variations studied here. Table |l] presents the 
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p,n,p,n 


STAR37 


422 ± 22 


719 ± 74 


1100 ± 110 


PHENIX38, 


453 ± 33 


674 ± 78 


954 ± 85 


L=0 


528 


897 


1310 


L=800 MeV/fm^ 


433 


714 


1027 


L=1.6 GeV/fm^ 


403 


652 


931 




406 


659 


945 


cl = 0.1 


433 


714 


1027 


cl = 0.2 


463 


772 


1116 


A-wr)/ s=Q 


408 


664 


957 




433 


714 


1027 


4nri/s—4 


449 


743 


1081 


Initially isotropic 


428 


695 


1012 


47r(C/s)max = 


462 


763 


1107 


47r(C/s)max = 2 


433 


714 


1027 


47r(C/s)max = 4 


418 


679 


983 


CGC IC 


447 


741 


1062 


Wounded Nucleon 


433 


714 


1027 


Collision Scaling 


482 


806 


1173 



TABLE I: The mean (pt) in MeV/c for central collisions for pions, kaons and protons. Only charged species were used 
in the PHENIX analysis, and only negative hadrons were used for STAR. Increasing the stiffness of the equation of state 
or the shear viscosity raises (pt) for heavier particles due to the corresponding increase in radial flow. The increase of 
radial flow with shear viscosity derives from the increase in Txx and Tyy relative to T^z at early times. By beginning the 
calculation with an isotropic stress energy tensor, radial flow is modestly reduced. Increasing the bulk viscosity lowers 
the effective pressure, and thus modestly reduces radial flow. Using coUisional scaling to set the initial energy density 
profile results in a more compact initial source, which then generates more radial flow. 

mean pt for pions, kaons and protons. Again, it should be emphasized that these calculations include all weak 
decays of hyperons and of the Kg- To some extent, these decay products are subtracted from experimental 
analyses, which might lead to the model predictions under predicting the mean pt for pions. However, the 
calculations also neglect symmetrization effects on the pion spectra, which should lower the pions mean pt- 
Unfortunately, the uncertainties in the experimental values in Table |T] are rather large, mainly due to the fact 
that experiments measure in a finite pt range. The experimental values for mean pt do agree with the default 
calculation, within the large experimental uncertainties. A more meaningful comparison, which would involve 
comparing actual spectra in the measured regions, is outside the scope of this study, but Table |l] is certainly 
sufficient for evaluating the sensitivity of the spectra to the various model parameters studied here. It should 
be emphasized that the sensitivity of spectra to the equation of state has been considered by numerous authors, 
and that in [35], sensitivity of the elliptic flow is also considered. 

A. Adjusting the Equation of State 

To study the sensitivity to the equation of state, we vary both the speed of sound and the width of the 
intermediate region, with the five equations of state being displayed in Fig. [5] For temperatures below 170 
MeV, or equivalently for energy densities below eh ~ 400 MeV/fm'^, the equation of state is that of a resonance 
gas. For the intermediate region, eh < e < £h + L, the equation of state has a constant speed of sound, 
P — Ph = Cg(e — eh)- Above the intermediate region, the speed of sound is set to — 0.3, to be consistent with 
lattice calculations at high temperature. 

The equation of state can be softened by either increasing L or decreasing the mixed-phase value of c^. Figure 
[6]displays source sizes for the three values of L: the default value of 800 MeV/fm'^, a soft value of 1.6 GeV/fm'^, 
and a hard value of zero. The default value was chosen to be crudely consistent with the behavior of lattice 
calculations which show a strong stiffening of the matter for energy densities rising from 1 — 1.5 GeV/fm^. The 
equation of state was also altered by adjusting the mixed- phase value of cl from the default value of 0.1 to 
either a stiffer value of 0.2 or to a softer value of zero, which would correspond to a first-order phase transition. 
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FIG. 5: (color online) Pressure vs. energy density for five equations of state: The default equation of state (solid line) 
assumed a constant speed of sound, — 0.1, for an intermediate range of energy densities, < e < + L, where 
L = 800 MeV/fm^. The speed of sound in the intermediate region was varied (dotted lines) to either 0.2 or zero. 
The latter choice corresponds to a first order phase transition. Keeping the default speed of sound, the width of the 
intermediate region was also varied (dashed lines) to either zero or 1.6 GeV/fm^. 

The femtoscopic effect of varying the speed of sound is shown in Fig. [7] 

As expected, softer equations of state lead to longer relative values of -Riong and -Rout relative to Rsidc- 
Whereas the increase in i?iong signals an increase in the mean emission time, the increase in i?out is indicative 
of a longer duration of the emission, or a more outside-in nature to the emission. The product of the three radii 
increases for softer equations of state. This is due to the increase in entropy associated with a softer equation 
of state (when compared at the same energy density) . Although these variations in the equation of state are 
rather strong, doubling i or c^, femtoscopic radii were affected on the level of 10%. Spectra are also affected 
by changes to the equation of state at the level of 10% as seen in Table [l] In particular softer equations of state 
lower collective radial flow which leads to lower values of the mean pt , especially for heavier particles. 



B. Adjusting Viscosities 



Even modest viscosities signficantly modify the stress energy tensor. It has been proposed that shear viscosity 
can not fall below the KSS limit, f] > s/Att [23]. According to Navier Stokes, at early times where the velocity 
gradient is 1/r, the KSS limit yields 

blTT \ 6Trl T J 

where the expressions involving P assumed a free gas equation of state, P = e/3. One expects for thermalization 
times near 1/2 fm/c where Tt ~ 1, that the correction to the longitudinal pressure is ^ 40%. If the viscosity is 
more than twice the KSS bound the value of Tzz can become negative. One expects a higher shear to accelerate 
the radial flow and result in lower values of i?iong and i?out/^sido [24^, as well as increased (pt) for heavier 
particles. These expectations have already been demonstrated by Romatschke [7 . 

The default calculation presented here assumes that the shear viscosity is twice the KSS bound. This would 
yield negative Tzz at early times, if not for the mapping described in Eq. (|6| which restricts the modification 
from shear to be less than the absolute value of the pressure. In the default calculation the initial condition for 
the stress energy tensor was set to this maximum value, with Tzz = and T^x ~ Tyy = e/2, consistent with 
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FIG. 6: (color online) Gaussian source dimensions for different equations of state. The equation of state incorporated a 
soft region where the speed of sound was set to = 0.1. The width of that region was L = 800 Mev/fm^ in the default 
calculation (circles) and is compared to L = (squares) and L = 1.6 GeV/fm^ (triangles). The softer equations of state 
yield larger values of i?out /i?iong . Experimental values are also depicted (stars). 



color-glass calculations [35] • For pure non-interacting classical fields, the anisotropy would be even larger as Tzz 
is negative, T^z — — e and T^x = Tyy = e. 

In addition to the default calculation, three modifications to the shear viscosity are illustrated here. First, the 
viscosity in the high density phase is set to zero. For the first modification the initial anisotropy to the stress 
energy tensor is also set to zero. For the second modification, the viscosity in the high-density phase is doubled 
relative to the default calculation to four times the KSS bound. Due to the ceiling imposed on the viscous 
modifications, the higher viscosity only matters for times greater than 1 fm/c. For the final modification, the 
default calculation is modified by setting the initial anisotropy of the stress energy tensor to zero. This mainly 
affects the expansion during the first one fm/c. Since the Israel-Stewart relaxation times tend to be ~ 1/2 fm/c, 
memory of the initial condition is lost after that point. 

The expectation for the femtoscopic radii are borne out by the results illustrated in Fig. [8] The default 
calculation differs from the zero- viscosity calculation by ^ 10%. In particular, the i?out/^sido ratio comes 
significantly closer to unity. It should be pointed out that the zero-viscosity calculation differs from many 
previous ideal hydrodynamic calculations in that the initial time was set to 0.1 fm/c, whereas several other 
calculations used either 0.6 or 1.0 fm/c, which would further increase the i?out/^sidc ratio. However, there is 
no physical justification for setting T^x — Tyy = at early times, thus such an initial state seems unwarranted. 
This is discussed in more detail in [4Dj. Modifying the initial anisotropy is similar in principal to altering the 
viscosity for early times. 

Bulk viscosity is only expected to be significant near the critical region. In particular, the condensed fields 
may not be able to keep pace with rapidly changing equilibrium values. This can lead to a peak in the bulk 



14 



.1.4 



^1.0 

o 



3 

o 

a: 6 



4 
8 



0) 



12 



cr 8 
4 



— H — £r ♦ ^ -c 



8 - 



0^ 6 - 



4 








100 200 300 400 500 
(MeV/c) 



FIG. 7: (color online) Gaussian source dimensions for different equations of state. Rather than varying the width of the 
soft region as was shown in Fig. k\ the speed of sound in the soft region was varied from the default value =0.1 
(circles) to = (triangles) and = 0.2 (squares). Experimental values are also depicted (stars). 



viscosity in the intermediate energy region [24j , which has been verified with analysis of lattice results |25j . 
The divergence of the velocity V • incorporates velocity gradients in all three directions is approximately one 
third (fm/c)~^ for r = 5 fm/c. For the default calculation, the peak value the ^ in the intermediate region 
is 2s/4tt. The magnitude of the effect is similar to what was derived in [24]. For such a velocity gradient the 
trace of the stress-energy tensor is modified by a substantial fraction. Doubling the bulk viscosity can make 
the Navier-Stokes value of (Tu) < 0. For the calculations performed here, the mapping procedure of Eq. ^ 
saturates the size of the change in {Ta) to be less than P. However, when combined with shear effects individual 
components can fall below zero. It should be emphasized that the non-equilibrium effects that generate bulk 
viscosity, mainly non-equilibrium fields, may be very poorly represented by viscous formalisms. First, in the 
transition region, responses may be highly non-linear, and secondly the field might not relax exponentially 
toward equilibrium as is assumed in Israel-Stewart treatments. Thus, the study here can really only point at 
the qualitative impact of bulk viscosity on the dynamics, and ultimately on the femtoscopy. 

The default calculation was modified twice, once by doubling the bulk viscosity and once be eliminating it. 
The difference of the three calculations, illustrated in Fig. |9] was modest. The effect on the stress energy 
tensor is similar to what one would get by softening the equation of state in the transition region. But unlike 
the changes to the equation of state investigated in the previous sub-section, these changes do not affect the 
pressure at either high or low density. Hence, the impact on observables tended to be modest. Increasing the 
bulk viscosity only affected the femtoscopy at the level of a few percent. The bulk viscosity helped amplify the 
magnitude of the pulse in the energy density created near the soft region. However, the pulse largely dissipated 
later in the collision. The bulk viscosity had a slightly more significant impact on the mean pt as seen in Table 
P^and in the initial energy density, which was adjusted to match dN/di], as seen in Table Ih] 
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FIG. 8: (color online) Gaussian source dimensions for three diflterent shear viscosities in the high-energy phase: Results 
from the default calculation with Anrj/s = 2 (circles) are compared to results using 4nrj/s — (squares) and Atttj/s = 4 
(triangles). Higher shear viscosities result in more rapid initial accelerations and smaller -Rout/-Rside ratios. Also shown 
is a modification of the default calculation where the stress-energy tensor was initialized as isotropic, rather than at the 
saturated values. Even though non-negligible variances resulted for the source radii, differences were found for the mean 
Pt seen in Table |l] Experimental values are also depicted (stars). 



Bulk viscosity had a visible impact on the smoothness of the energy density profiles. Larger bulk viscosities 
appeared to lead to jagged and unstable profiles in the intermediate region. We speculate that this is driven 
by the fact that bulk viscosity effectively pushes the pressure vs. energy density to behave non-monotonically, 
which could give regions where d{Tn)/de, which is the effective speed of sound squared, is zero or negative. It 
would be interesting to know whether such instabilities would appear in a more physically grounded treatment 
of non-equilibrium effects, such as one where the dynamics of non-equilibrium fields were treated in parallel to 
the hydrodynamic treatment by solving a coupled Klein-Gordon equation. 

The significant sensitivity of the final femtoscopic source sizes to acceleration during the first one or two 
fm/c might seem surprising. The importance of early acceleration can be likened to an Olympic sprint, where a 
head start of a few tenths of a second results in a difference of several meters at the end of the race. For this 
reason, it is imperative to understand the bulk properties, e.g., the stress-energy tensor, of matter even before 
thcrmalization . 



C. Adjusting Initial Conditions 



For a rotationally and boost invariant calculation the choice of initial conditions involves choosing the initial 
transverse density, the initial stress energy tensor, and the initial time. For all calculations presented here, the 
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FIG. 9: (color online) Gaussian source dimensions for three different bulk viscosities. Bulk viscosities were set to zero 
outside the intermediate-energy density region near Tc, and varied linearly from zero to a maximum value in the center 
of the region. The three values, 47r(^niax/s= (squares), the default value 2 (circles) and 4 (triangles), resulted in nearly 
identical radii. Somewhat larger differences were observed in the mean pt values from Table |l] and in the initial energy 
densities from Table [ll] Experimental values are also depicted (stars). 



initial time was chosen to be 0.1 fm/c. Since the development of collective flow at such early times is driven by 
the uiitial stress energy, there is no reason to pick a later time, unless there were reason to expect Tj-r^ and Tyy 
to be zero at early times. Since a little more than 0.1 fm/c is required for the nuclei to pass one another, and 
given that this is already an extremely short time relative to the overall expansion time, there is no motivation 
to pick an earlier time. Variations of the initial stress energy tensor, and in particular variations to the initial 
anisotropy of Tij were considered in the previous subsection along with variations in the viscosity. 

Three variations of the initial energy density profile were explored. The default calculation was that of the 
wounded nucleon model |41| . In this calculation the probability of a nucleon interacting is calculated as unity 
minus the probability it survives without interaction. The survival probability is calculated assuming the particle 
travels through the Woods-Saxon nuclear profile with a 40 mb cross section. For the thick part of the nucleus, 
this approaches participant number scaling. An alternative is collision scaling, where the energy density at a 
transverse coordinate {x, y) is proportional to TaTh, where the thickness function T is calculated by integrating 
the density of a nucleus over the z coordinate. The third profile explored here is for the color-glass profile used 
in "42]. In that case the energy density is chosen proportional to the minimum of Ta and T{,. For all three 
profiles, the thickness functions were found by convoluting two Woods-Saxon profile whose centers differed by 
an impact parameter of 2.21 fm, corresponding to the 5% most central collisions of Au-|-Au. The density profile 
was then averaged over the azimuthal angle to generate an approximate radial profile. For every calculation, 
the profiles were renormalized so that the resulting dNch/drj was 691 ± 5, consistent with [55] . 

The collision-scaling profile was the most compact of the three attempted here, and resulted in the largest 
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FIG. 10: (color online) Gaussian source dimensions for three different initial energy density profiles: the wounded 
nucleon model (circles) is used for the default calculation, and results in a less compact source than either the color-glass 
inspired model of [35] or collision scaling. More compact sources are more explosive and lead to lower Rout/Rside ratios. 
Experimental values are also depicted (stars). 



radial flow. This profile resulted in higher (pt) for protons and lower values of -Rout/^sido as seen in Fig. 10 
The least compact profile was the default calculation which was based on the wounded nucleon model. 

The initial energy density at the extreme periphery of the collision should be driven by coUisional scaling 
since one can ignore the possibility of multiple collisions in that limit. None of the three profiles employed here 
obey this constraint, as the overall normalization is scaled in order to match the experimental value dNch_/drj. 
If this constaint were replaced by the constraint that the density profile behaved correctly for small Ta and T^, 
the experimental dNch/drj could instead additionally constrain the remaining parameters. The experimental 
dNch/drj has been used to argue that the profile appears more driven by participant scaling than collision scaling 
[35j . However, these analyses are exceedingly simple and neglect many aspect of the expansion such as viscosity 
or longitudinal work. Since the temperature is not constant across the profile, the relation between entropy 
and energy densities is not linear, thus different profiles are reached if one believes the energy density should 
follow a given scaling vs. the entropy density, which is more related to the multiplicity. This underscores the 
importance of performing a global analysis of all observables, including elliptic flow, which is also affected by 
the choice of proflle shape [42] . 

It is unfortunate that the initial energy density is not itself an observable. After adjusting the initial energy 
density for each parameter set to match the experimental dNch/drj, the final state observables tend to change 
at the 10% level or less for all variations studied here. However, the initial energy densities vary by more than 
a factor of two for these calculations as seen in Table Ull This variation has little to do with the variation of the 
asymptotic transverse energy for fixed multiplicity, which is also a 10% effect. The variation largely derives 
from differences in the longitudinal work in the expansion. The work is proportional to both T^z and the time 
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width of soft region in EoS 


eo (GeV/fm^*) 


L=0 


150 


L=800 MeV* 


114.5 


L=1.6 GcV 


104.5 


stiffness of soft region in EoS 




cl = Q 


107 


cl = 0.1* 


114.5 


d = 0.2 


124.5 


shear viscosity in parton phase 




Itt 77/5=0 


289 


Itt 77/5=2* 


114.5 


Itt 77/5=4 


106.5 


initially isotropic init. cond. 


148 


max. bulk viscosity in soft region 




47r(;^max/s = 


124 


47rCmax/s = 2* 


114.5 


47rCmax/s=4 


109 


initial density profile 




CGC IC 


136 


Wounded Nuclcon* 


114.5 


Collision Scaling 


180 



TABLE II; Initial central energy density at tq = 0.1 fm/c, in GeV/fm^. Values, which were adjusted to match 
experimental values of dNch/drj, vary by more than a factor of two, largely due to differences in longitudinal work. The 
default calculation (noted by *) is varied in five ways. 

over which the system expands. For large viscosities, softer equations of state, or for initial conditions where 
Tzz is small, the longitudinal work is reduced, thus requiring smaller initial energy densities to attain a given 
final condition. For shear anisotropics T^x and Tyy are also enhanced, which accelerates the expansion. This 
results in a reduced time for the reaction, which also reduces the longitudinal work. Changing the shape of the 
initial profile also changes the energy density, mainly for the trivial reason that a more compact initial profile 
requires a higher energy density to produce the same dNch/drj. 

IV. SUMMARY AND OUTLOOK 

The principal conclusion from these investigations is that the femtoscopic data from RHIC can be repro- 
duced to within 10% with models combining viscous hydrodynamics and hadronic cascades. In particular, the 
i?out/-Rsidc ratio can be brought down close to unity. The failure of previous models appears to derive mainly 
from three shortcomings, all of which are related to under-predicting the explosivity of the collision. First, the 
equations of state were often too soft, using a first-order phase transition. A stifFer equation of state is more 
explosive, and can lower the i?out/^sido ratio. Secondly, previous treatments ignored acceleration before the 
thermalization time. From general arguments involving conservation of energy and momentum in the equations 
of motion of the stress-energy tensor, it should be clear that thermalization is not required for acceleration. 
In fact, longitudinal classical fields, which are far from equilibrium by definition, result in strong transverse 
acceleration. Finally, the previous treatments were based on ideal hydrodynamics. The effects of shear, as 
already demonstrated in [7] , increase the transverse pressure relative to the longitudinal pressure at early times, 
which of all the variations considered here, appears to be the most important. Bulk effects, were manifest in 
the final mean pt, but made remarkably little difference in femtoscopic radii. Previous treatments overpredicted 
the i?out/^^side ratio by 40% or more |30j . a result confirmed if we run this model with a softer equation of state, 
without viscosity, and delaying transverse acceleration for the first fm/c. 

It would appear that improving models in all three areas mentioned above is required for rectifying the short- 
comings. The default calculation, which includes all three such effects, provides a reasonable description to the 
data. Without including longitudinal acceleration, which requires a three-dimensional model, it is unreasonable 
to expect better agreement and it is probably not meaningful to try to better fit the data by adjusting parame- 
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ters. An additional area of uncertainty documented here conies from the choice of the initial profile, as a more 
compact source results in a more explosive source. One could reduce any of the three effects mentioned in the 
previous paragraph, then compensate for them by adjusting the initial density profile. 

A second impression generated by this investigation is that it appears impossible to disentangle the various 
uncertainties mentioned above by focusing only on two-pion intcrferometry. Spectra are sensitive to the same 
model features studied here, as evidenced by the mean pt values listed in Table |l] Elliptic flow observables, 
which require a higher-dimensionality model than used here, can also be used to assist in understanding the 
collision. Hopefully, different observables will be relatively more sensitive to different facets of the model. Then 
by performing a coordinated analysis of numerous classes of observables, one should be better able to answer 
specific question and determine specific parameters. These analyses should also incorporate a greater set of 
femtoscopic measurements, which we list here: 

• Femtoscopy using particles other than pions. Heavier particles are more sensitive to collective flow due 
to their lower thermal velocity. Correlations between a heavy particle and a light particle, e.g. pn, are 
especially sensitive to the patterns of collective flow. 

• Non-central collisions and collisions at different energies. Measurements have already been made as a 
function of the direction of the pair's momentum relative to the reaction plane |28j . This information is 
rich in detail, but the meaning of the information is not yet understood. Additionally, there exists data 
at different beam energies. By studying the response to changes in the initial energy density, without 
changing the size, one should gain some leverage for disentangling some of the issues mentioned above. 

Finally, it should be emphasized that, although the large discrepancy with the i?out/^sidG ratio has been 
eliminated, none of the variations studied here provided a completely satisfactory reproduction of the pt de- 
pendence of source dimensions. The data showed a modest fall of the ratio for higher pt, which combined with 
the constraint that the ratio must be unity for pt = 0, gives a non-monotonic behavior. Although the rise and 
fall are only of the order of 10%, the model calculations all showed monotonic behavior with pt- Furthermore, 
the model calculations tend to over-predict i?iong at low pt. This might be due to the assumption of boost 
invariance, which if relaxed, should provide corrections in the direction of the data. Finally, conspicuous by 
its absence, has been a comparison of the A factors, which represent the fraction of pairs that are correlated. 
The model calculations over-predicted these factors, but without a better understanding of experimental details 
about acceptance of weak decay products and particle mis-identification fractions, one can not as yet draw any 
conclusions. 
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